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Abstract 

A globally converging numerical method to solve coupled sets of non-linear integral 
equations is presented. Such systems occur e.g. in the study of Dyson-Schwinger 
equations of Yang-Mills theory and QCD. The method is based on the knowledge 
of the qualitative properties of the solution functions in the far infrared and ultra- 
violet. Using this input, the full solutions are constructed using a globally conver- 
gent modified Newton iteration. Two different systems will be treated as examples: 
The Dyson-Schwinger equations of 3-dimensional Yang-Mills-Higgs theory provide 
a system of finite integrals, while those of 4-dimensional Yang-Mills theory at high 
temperatures are only finite after renormalization. 
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1 Introduction 



Integral equations appear in many areas of physics, and they are ubiquitous 
in non-perturbative field theory. Most prominent examples are the Dyson- 
Schwinger equations (DSEs). These equations form coupled sets of non-linear 
integral equations, which must be solved self-consistently to make progress 
beyond resummed perturbation theory. Solutions to these equations thus lie at 
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the heart of many truly non-perturbative calculations. E.g. in the low-energy 
regime of QCD, the study of the DSEs has provided a wealth of physical 
insight, see e.g. ref. [1] for a review. To solve these equations, fast and reliable 
algorithms are desirable, which minimize the required prior knowledge. 

While the theory of single linear integral equations with well-behaved integral 
kernels is well understood, see e.g. [2], this is not the case for typical equa- 
tions in field theories. These are normally coupled sets of non-linear equations 
with integrable, but singular integral kernels and not necessarily well-behaved 
solutions. 

In many problems where solutions are only searched on a compact manifold, 
discretization of the integral equations allow to solve them by fix-point iter- 
ation methods. This is especially the case in non-relativistic settings, such as 
condensed matter problems. In fully relativistic settings, such solutions are 
plagued by finite volume effects and discretization errors for extremely large 
and extremely small momenta. This is in some cases of no importance, but in 
the case of non-abelian gauge theories, these regions are of special interest, as 
they are the origins of non-perturbative phenomena and of matching to per- 
turbation theory. Nonetheless, methods using discretization provide a wealth 
of valuable informations also in these cases, see e.g. refs. [3,4,5,6]. 

The aim here is to construct a globally convergent method applicable to such 
problems. In this paper section 2 will layout the DSEs to be solved and thus set 
the mathematical frame 1 . The numerical method, including the implementa- 
tion, will be discussed for two different systems. A finite system is obtained for 
3-dimensional Yang-Mills-Higgs theory. It is treated in section 3. In section 3 
also an explicit example will be shown. A renormalized system is obtained for 
finite-temperature 4-dimensional Yang-Mills theory to be discussed in section 
4. In section 5 some concluding remarks are made. Some technical details are 
deferred to an appendix. A thorough discussion of the physics case can be 
found in refs. [7,8,9]. 



2 Dyson-Schwinger equations 

The complete content of a theory is given by the corresponding Green's func- 
tions. The Dyson-Schwinger equations [10] determine the Green's functions 
of a theory through an infinite set of coupled, non-linear integral equations. 



This section includes a short description of the derivation of the equations to be 
solved. If the reader is exclusively interested in the numerical method, this section 
can be mostly skipped, except for the systems of equations to be solved. These are 
given in (9-11) for the finite system and (15-17) for the renormalized system. 



2 



A theory is fully described by the Lagrangian density £(0 a ,<9 M a , ...) of the 
fields <p a and their derivatives, where a is a generic multi-index. The Green's 
functions 2 are then given by functional derivatives w.r.t. the fields <p b of the 
identity 



Here d is the number of space-time dimensions. See refs. [1,11] for a detailed 
introduction to DSEs. 

It is in general not possible to solve all DSEs simultaneously, and it is therefore 
necessary to truncate the system to a finite number of equations. This entails 
nearly always violations of internal symmetries such as gauge symmetries, and 
it is an essentially non-trivial task to compensate or estimate these effects. 
For a detailed study of this problem, see e.g. refs. [1,6,7]. The consequences of 
these truncation and how to control them is a physical issue, which will not be 
treated here. The implementation of appropriate measures to deal with this 
problem can have impact on the numerical method. Furthermore, in quantum 
field theory, terms in the DSEs usually diverge, and must be made finite 
by renormalization [12]. Although this is again a physical issue, it has also 
consequences for the numerical method. Therefore some numerical aspects of 
both problems will be discussed briefly in section 4. 

The presented numerical method has been applied successfully to two systems. 
One is 3-dimensional Yang-Mills theory coupled to an adjoint Higgs as the 
infinite-temperature limit of 4-dimensional Yang-Mills theory [7,8]. The second 
is 4-dimensional Yang-Mills theory at high temperature [7,9]. The first system 
is finite and the numerical problem thus will be treated first in section 3. The 
second requires renormalization and the necessary alterations will be discussed 
in section 4. In the following subsections, the Dyson-Schwinger equations of 
each system will be shortly introduced. 

2.1 3-dimensional Yang-Mills theory coupled to an adjoint, massive Higgs 

The infinite-temperature limit of 4-dimensional Yang-Mills theory is described 
by a 3-dimensional Yang-Mills theory coupled to an adjoint, massive Higgs 
field [13]. For physical reasons the choice of Landau gauge is advantageous, 
which will be made here exclusively. The Lagrangian of this system in Eu- 
clidean space is then given by [8,13] 



2 The discussion here is restricted to an Euclidean space-time, but can directly be 
transferred to any other metric. 




(1) 



£ = \f^f; v + c a d,Dfc b + l -{Df<t> b D a ;4> c + ™ 2 h <t> a <t> a ) + (2) 
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with the field strength tensor F* v and the covariant derivative D ab defined as 

C = d»K - d u A; - g^ ahc A\Al (3) 
Dt = 8 a % + g 3 f abc A^ (4) 

A ^ is the 3-dimensional gluon field, c a and c a are the Faddeev- Popov ghost and 
anti-ghost fields, (fi a is the Higgs field, g 3 is the dimensionful coupling, which 
is 1 in appropriate units, and m h ~ g\ is fixed by Monte Carlo calculations 
on discrete space-time lattices to be 0.88gf [14]. The constants f abc are the 
structure constants of the gauge group. In general, contributions stemming 
from the Higgs field could be divergent. However, this is not permissible when 
viewing the Lagrangian (2) as the infinite-temperature limit of 4-dimensional 
Yang-Mills theory. Therefore leading order perturbation theory requires that 
h is fixed to 

where Ca = fabcf abc is the adjoint Casimir of the gauge group. 

The simplest non-vanishing Green's functions in this theory are obtained by 
deriving the identity (1) once more with respect to the fields. The inverse of 
these Green's functions are called propagators, and there is one for each field. 
These propagators D G , D^, and D H of the ghost, the gluon, and the Higgs, 
respectively, can be parameterized by three scalar functions, so-called dressing 
functions, G, Z, and H as 



W) = =f2 (6) 
D M =( S „-?f)^ (7, 

iW) = ^, (8) 

respectively. The dressing functions G, Z, and H have to be positive semi- 
definite by the Gribov condition [8,15]. Using these dimensionless quantities 
instead of the propagators improves the numerical stability for large p sig- 
nificantly, as otherwise trivial kinematic effects would require an enormous 
precision. 

The truncated DSEs are then obtained by standard techniques [1] in a straight- 
forward but tedious way [7]. The tensor equation for the gluon propagator 
is contracted with an appropriate projector to yield a scalar equation. The 
projector is parameterized by a variable (, on which therefore the integral 
kernels depend. This yields the DSEs for the ghost dressing function G, the 
gluon dressing function Z, and the Higgs dressing function H 
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= 1 + |^/ d0dqA T (p, q)G{q)Z{p - q) - J-j (9) 



= 1 + + T HG + T HH + 



p2 . - . - . ( ^ / d9dq^N 1 (p : q)H(q)Z(p + q) 
+N 2 (p,q)H(p + q)Z(q^- 1 ±- ) (10) 
= 1 + T G " + T GG + |^ | (?)(?((?)(?(? + q) 

+M L (p, q)H(q)H(p + q) + M T (p, q)Z(q)Z(p + g)) - ^-y, (11) 

where is the angle between p and g. The tadpoles T lJ are not independent 
functions but compensate any divergencies and thus render the equations (9- 
11) finite [7]. The integral kernels 3 A T , N u N 2 , R, M T , and M L as well as 
the tadpoles are listed in appendix A. Trivial factors, such as the integral 
measure, have been included into the kernels. The kernels all have a very 
similar structure, being rational functions of the momenta. All these kernels 
furthermore contain integrable singularities of the type \/{p — q) m with some 
positive m. 



2.2 Finite temperature Yang-Mills theory 



The second example is finite temperature 4-dimensional Yang-Mills theory. In 
the Matsubara formalism [16] it is described by the (Euclidean) Lagrangian 

c = \f; v f; v + c a d,Dfc\ (12) 

which is very similar to (2). The only difference to the previous example is re- 
placing (yf 3 by the 4-dimensional dimensionless coupling g± in the field strength 
tensor (3) and the covariant derivative (4) and dropping the Higgs field. Using 
the Matsubara formalism [16] the DSEs of the vacuum [1] can be extended 
to finite temperature [7]. The propagators are again the most simple Green's 
functions and can be parameterized by also three functions 

Da® = ^#^> (13) 

P A 



3 The kernel Mt in addition depends on a parameter Ss g . It is fixed throughout to 
1/4 and parameterizes the constructed three-gluon vertex. Furthermore this induces 
an additional non-trivial dependence of Mt on G and Z, see appendix A and ref. 
[8]. 
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for the ghost propagator and similarly for the gluon propagator as [16] 

DfM = PrM^l + P^M^P-. (14) 

where and Pl are projectors transverse and longitudinal w.r.t. the heat- 
bath [16]. The coincidence of the names of the dressing functions to those of 
the previous example is intentional. There is an exact correspondence in the 
limit of infinite temperature. 

The truncated DSEs for the dressing functions G, Z, and H are obtained in 
the same way as before and read 



Q-Z + 9iT ° A 

3 w 



00 r / 

£ J d9dq[A f T (p,q)G(q)Z(p-q) 



2 

n=—oo 



+A[(p,q)G(q)H(p-q))-— (15) 
= ZZ 3L + Tf HG + Tf HH + 9 ^ ± [d6dq(pf(p,q)G(q)G(p + q) 

( 27r ) n=~oc J V 

+7V/(p, q)Z(q)Z(p + g) + N{(p, q)H(q)Z(p + g) 

+JV/(p, (p + q)Z(q) + N*(p, q)H(q)H(p + g)) - ^ (16) 

(ATTj n= _ 00 J 



(27T) 

+M/(p, q)H{q)H{p + q) + M((p, q)H{q)Z{p + g) 
+M/(p, 9 )tf(p + g)Z(g) + M£(p, q)Z(q)Z(p + q) 



+ ^^~ r L " J " ^y- (17) 

The temperature is denoted by T. The summation is performed over all in- 
tegers n. In the Matsubara formalism the component go is discrete and given 
by g = 2nTn, called the Matsubara frequency. Also the external p is thus 
discrete. 

In the present case, the sum is truncated, and only the range [-N +1,N — 1] 
is included. Close inspection of the equations reveals further that the dressing 
functions can only depend on \p \ and \p\. The corresponding symmetry, under 
Po — > — Po, is used to reduce the number of equations significantly. Therefore 
3N independent functions have to be determined. The largest system treated 
numerically so far is iV = 21. 

The kernels A f T , A f L , R f , Af{,, Af(, M$, M[, P f , N{., N{, N$, and N[ are 
rather lengthy, and will not be displayed here. They can be found in ref. [7]. 
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The kernels are very similar in structure to their 3-dimensional pendants listed 
in appendix A. The same applies to the tadpoles T-™, with the notable ex- 
ception of th T^ H K These are much more similar to the transverse expressions 
T fGj then in the case of 3 dimensions. Especially these are now also deter- 
mined by integrals, see ref. [7] again for details. The finite-temperature theory 
is renormalizable. Thus explicit wave-function renormalization constants Z 3 , 
Z 3 l, and Z 3 t have been introduced 4 . 

The variables ( and £ are introduced by contraction of the equation of the 
gluon propagator with two different projectors [7,9] as before. These variables 
have been introduced to study the consequences of the truncation. For po = 0, 
equation (16) is only superficially dependent on £: since all integral kernels 
are proportional to £, the dependence can be divided out for £ ^ 0. 

The renormalization constants are defined at a subtraction point s, which can 
be chosen arbitrarily, but for numerical reasons should be different from 0. 
The explicit implementation of the renormalization prescription [7,9] for the 
DSE of a dressing function F = G, Z, H with self-energy contributions / 

' 1 + /(P) (18) 



F{p) 

is then for the corresponding wave function renormalization constant Z 3 

1 + 5Z 3 + I(p) (19) 



F(p) 

5Z 3 = -I(s) (20) 
Z 3 = 1 + 5Z 3 (21) 

and in the case of H(0, \p\) 



l + 5Z 3L + °— + I(p) (22) 



H(0M p 2 

Sm 2 =m 2 - lim p 2 1(p) (23) 

, Z „ = _ /W + «2«=*££M (24) 

where m r = m 3( i is the renormalized mass, and 5 is an infinitesimal parameter, 
different from zero only for numerical reasons. It is taken to be the same value 
as the infrared integral cutoff introduced in the next section. The difference 
compared to performing the subtraction at is negligible. Note that due to the 
truncation of the Matsubara sum the wave function renormalization constants 



No vertex renormalization occurs within this truncated system. 
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depend on p . This dependence vanishes for N — > oo [9], but has to be taken 
into account in the numerical treatment presented here. 



3 Numerical method 

The numerical technique to solve the system of equations (9-11) and (15-17) 
are very similar. Renormalization in the 4-dimensional system (15-17) requires 
some extensions. Therefore here only the finite system will be treated. The 
modifications for the renormalized system will be discussed in section 4. 

The earliest approach to solve the 4-dimensional equivalent (without the ad- 
ditional Higgs field) of system (9-11) was based on a direct iteration method 
[17]. This has afterwards been replaced by a method based on a local New- 
ton iteration [6,18,19]. The latter approach used as starting inputs for a local 
iteration results obtained from a discretized version of the problem [3]. The 
method presented here generalizes this ansatz to provide global convergence 
and remove the necessity of input data. For the sake of completeness the full 
method will be described here, although several elements have been taken 
from the previous one. 

The basis for the numerical treatment is an analytical investigation of the 
asymptotic behavior of the solutions. This will be discussed in subsection 
3.1. An important ingredient is the representation of the solutions using their 
asymptotics. This and the integration method will be discussed in subsec- 
tion 3.2. The heart of the method is the iteration procedure, which will be 
described in subsection 3.3. Some notes on the efficient implementation are 
given in subsection 3.4. An explicit calculation with a typical set of numerical 
parameters is given as an example in section 3.5. 

3.1 Analytical foundation 

It is important to gather as many analytical properties as possible to obtain an 
appropriate representation of the solutions. Especially using as few parameters 
as possible permits to obtain maximal numerical efficiency 

The searched for positive semi-definite solution functions are continuously dif- 
ferentiable and analytical in the momentum range (0, oo). For physical reasons 
no poles at intermediate momenta are expected. However the solutions are not 
bounded from above. Thus divergences at p = are possible. Therefore the 
infrared asymptotic behavior is of special interest. Favorable, it turns out that 
the leading asymptotic behavior can be obtained analytically [7,8,20]. 
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In the infrared, power-like ansatze can be used [8,20] 5 , 

limF(p) = A F (p 2 )~ eF . (25) 
p— >o 

This is possible due to the singularity structure — q) m of the integral 
kernels, which lead to a domination of the integrals by momenta of order q ps p. 
It is hence allowed to replace the dressing functions by the ansatze (25). It 
is then possible to solve the massless integrals analytically [8,20]. The Higgs 
propagator is dominated by the renormalized mass term in equation (10) and 
thus has the unique exponent e# = —1- This also determines Ah uniquely [8]. 
There are two possible solutions for the ghost exponents. One of the solution 
is e G = 1/2, while the other depends on ( and varies as e G e (1/4, 1/2]. The 
gluon exponent ez is not independent, but related to e G by [8,20] 

ec= 4( ez+ £>- (26) 

The coefficients Aq and Az cannot be determined uniquely. However, they are 
related by [8] 



1 _ C A gl 2 4 ^-!)(2 + 2e G (C - 1) - C)r(2 - 2e G ) sin 2 ^ 



A 2 G A Z (4vr)§ cos(27re G )(e G -l)e G r(|-2e G ) 
1 

= Ae G ,CY 



Furthermore, in the ultraviolet, p 3> gf , an analytic solution can be obtained 
by perturbation theory due to asymptotic freedom [8]. The dressing functions 
are then given by 

hm F(p) = 1 + ^M^, (28) 

The wf are the positive constants wq = 1/16, wz = 9/64 and wh = 1/4. 
Note that in contrast to 4 dimensions these are already the resummed solu- 
tions, as leading order and leading order resummed perturbation theory in 
3-dimensional Yang-Mills theory are the same [8]. This concludes the list of 
known analytical properties of the solutions. 



3.2 Expansion and integration 



The first step is to describe the solutions in terms of known functions. The 
solutions are due to their ultraviolet behavior neither integrable nor square- 
integrable. This prohibits an expansion in any orthogonal set of polynomials. 

5 Note that this is very similar to the case of 4 dimensions [21]. 
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This problem is reinforced by the infrared divergence of the ghost dressing 
function. 

Therefore originally [18,22] the dressing functions were expanded as 

F(p)=e(e-p)D I F (p) 

+9(p - e)0(A - p)D%(p) + 6(p - A)D u F {p) (29) 

N c h 

\n{p?{p))=Y,<%T i {M{p)) (30) 

3=0 

where Dp denotes the analytical infrared solution (25) and Dp the analytical 
ultraviolet solution (28). The index % runs through 1, 2, 3 corresponding to 
G, Z, and H, respectively. The range [e, A] is covered by a numerical expan- 
sion 6 in N c h Chebychef polynomials Tj. Good results can be achieved with 
— 30 — 40. For special purposes like computation of Schwinger functions 
or thermodynamic quantities [8] N c h = 100 — 500 are used. 

The typical expansion range in the case of 3 dimensions is chosen to be [e, A] = 
[10 _2 gf , 2- 10 3 (7jj]. The function M in (30) is a suitable mapping function which 
maps the momenta to the domain (—1, 1) of the Chebychef polynomials. Since 
the dressing functions F = G, Z , H vary on logarithmic scales, the mapping 



function is chosen accordingly as 

M(p) = A + B ln(p) (31) 
where A and B are chosen such as to yield 

M- 1 (e) = z (32) 
M~ 1 (A) = z Nch , (33) 



where M 1 is the inverse function of (31) and Zi are the ith zero of the Cheby- 
chef polynom of order N c h- 

To increase the numerical stability, (29) is altered to 



6 Expanding the logarithm instead of the function itself leads to a considerable 
smoothing at large momenta and the number of needed Chebychef polynomials 
drastically decreases. 

7 An alternative is a conformal mapping like z/(z + l) which removes the necessity 
of a numerical upper and lower cutoff of the integrals. However, such a mapping 
does not provide any advantage beyond formal elegance and in addition did not 
provide as well spaced evaluation points as the logarithmic mapping. 
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V 
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V 




Higgs (H) 


l 








Table 1 

The coefficients for the fitting functions (35-37) in terms of the scale of the problem 
v and of the renormalized Higgs mass m r [8] . 



D F (p) = 6(e-p)D F (p) 



+ 0(p - e)0(A - p)Dp (p)f F (p) + 9(p - A)D u F (p) 



(34) 



where the fp interpolate between the qualitative behavior of D F and D F . 
Thus for divergent quantities like the ghost dressing function 



faip) = 1 + 



a G b G +p- 2e Q 

p~ 2e G Cq + p 



(35) 



with suitably chosen fitting constants aa, be, cq- This especially improved the 
quality at large p, as it is notoriously difficult to fit with a polynomial onto a 
constant. Similar for converging quantities as the gluon and the Higgs the fits 



&hP 2eH 

fr<P> V- + H (37) 

are used, respectively Note that it turns out to be more stable to neglect 
in the case of the Higgs the subleading contributions in (28) fitted by c#. 
The coefficients are chosen to resemble the typical scale v for the system 
investigated. In case of the 3-dimensional theory this is v — g\ while in case of 
the 4-dimensional system v = g\T . The coefficients can then be read off table 
1 for the 3-dimensional theory. In the 4-dimensional theory a more simplified 
version of (35-37) is used, omitting the perturbative tail. For the hard modes, 
Po 7^ 0, essentially a fit of type (37) is used. 

Additionally it is necessary to choose an integrator. For the radial integral, 
after applying the logarithmic map (31) to the integration variable, Gauss- 
Legendre integration is sufficient [22]. Since the integrals are finite, a lower 
and upper numerical integral cutoff has been used, and the integration has 
been performed in [5,0], with 5 e and A f2, typically [10~ 5 gf, 10 7 p|]. 
The remaining problem to be dealt with is the singularity structure of the 
integral kernels and of the ghost dressing function G. These appear at \q\ = 
and p = ±q, depending on the equation. The \q\ = singularity is cured by 
imposing the lower cutoff. Here the logarithmic mapping of the integration 
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variable provides a sufficient sampling of the integrable singularities arising. 
More problematic is the singularity due to the coincidence with the external 
momentum. It can be removed by splitting the integral in two parts, up to 
and above p. Again the logarithmic spreading provides sufficient sampling of 
the integrable singularity. As p will only be chosen from [e, A] the integral is 
finally performed as 

oo e p A Q 

hhhhl- < 38 > 

5 e p A 

The angular integral is straightforward and can be done by normal Gauss- 
Legendre integration. However, it is necessary to perform the angular integral 
before the radial one, as some divergences cancel only due to the angular 
integration. These cancellations also require a sufficient sampling, leading to 
N a = 80 points in angular direction and N s + N £ + N p + N A = 60 + 100 + 100 + 
60 = 320 points in radial direction as typical values. 

3.3 Micro-, Macro, and Super-cycles 

Originally, solutions for the set of coefficients c^- were obtained using a local 
Newton iteration procedure [18]. The presented method here uses instead a 
global Newton method with back-tracking [23]. This approach is very insensi- 
tive to the starting value, and a typical starting guess is c^- = 0. The necessary 
3N c h functions to determine the coefficients are obtained by evaluating the 
system of equations (9-11) at the zeros of T^ ch . At these points the expansion 
in Chebychefs is exact. This also entails that all evaluation points stem from 
[e, A]. This makes it even more necessary to know the behavior of the functions 
outside this domain analytically to perform the integrals in the intervals [5, e] 
and [A, Q]. 

Therefore the numerical problem was to solve the 3N c h non-linear equations 
for the coefficients c^- and determine the coefficient Aq. These equations can 

— * 

be collected generically into an equation vector E, with components defined 

as 

E {i - 1)Nch+j = r{ Pj ) + £/(F i ,F I ,p i ) - (39) 

1=1 *i\Pj) 

with Fi = G, F 2 = Z, and F 3 = H and i — 1,2, 3. The index j runs from 1 to 
N ch , I are the integrals, and r the remaining tree- level terms. Of course, an 

— * — * 

exact solution of the system satisfies E — 0, and thus the aim is to find such 
a solution. 

New values c'^ for the coefficients are then obtained iteratively by [23] 

^ = 0- \.r l c, (40) 
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where c is the vector of coefficients Cy with the components (c)(i-i)N ch +j — Qj- 
J is the Jacobian defined as 

^ - ^ <«> 

The stepwidth is determined by the backtracking parameter A. It is chosen 
such as to minimize \E\, where a simple quadratic interpolation has been used 
[23]: A first trial step is a full local Newton step (A = 1). This is accepted, if 

— * 

it decreases the previous value \E\ = r . If this is the case, the algorithm has 
already moved so close to the correct solution that the faster local algorithm 
can be used (which provides quadratic instead of the linear convergence of the 
global version [23]). 

If this not leads to a decrease, a trial step of size o\ < 1 is performed. If 
this again provides no reduction of \E\, the quadratic interpolation is used: 
Given the step-size of the current trial step A 2 and of the previous trial step 
Ai, leading to values r 2 and r 1 of \E\, respectively, a new trial value of A is 
determined by a quadratic interpolation using the first two terms of a Taylor 
expansion of \E\ to construct [23] 



Pi = 



ro(A| 



A?) - \jr 2 + Aln 



AiA 2 



2 ~ Ay'A L > 



P2 



(2AiA 2 )(A 2 (ri - r ) - Ai(r 2 - r )) 
Ai — A 2 



Thus if P2 is greater than zero a descent direction is found, and the new value 
A = —pi/p2 minimizes the Taylor series truncated at the quadratic order. 
However, the quadratic approximation may not be sufficient, leading to a 
step-width outside the interval (0, A 2 ] or may be too close to or A 2 to be 
reliable. Thus if the new A lies outside the range [<7oA 2 , <7iA 2 ], with 00 < <7i, 
the new value of A is taken to be the corresponding border value, e. g. if the 
new value of A would be larger than <7iA 2 it is reduced to OiA 2 . If p 2 is positive 
or zero, a quadratic approximation is not reliable, and the new trial step is 
selected as A = (XiA 2 .With this new value of A, a new trial step is performed. 

— * 

If still no decrease is found, this procedure is repeated until \E\ is decreased, 
and the trial step is accepted and (40) is performed. This is repeated until 
convergence is achived, which will be discussed below. 

The method fails to converge if A becomes less than a predetermined value 
oiq. In the present case «g is essentially determined by empirical knowledge, 
but in general cases should be at least of the order of the machine precision. 
The values of a and cr 1 are tuning parameters for the algorithm, a specific 
choice is presented in subsection 3.5. 
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A significant problem 8 is the determination of the coefficient Aq of the in- 
frared asymptotic behavior (25) of the ghost dressing function G. Trying to 
incorporate it into the global Newton method failed. Hence this part is moved 
to an outside iteration: after each global Newton iteration, a new value of Aq 
is determined by requiring that the infrafred solution D G connects continu- 
ously to the intermediate function JqDq. In general this leads to a continuous 
but not to a continuously differentiable function until convergence has been 
achieved. 

Therefore A G has not yet its final value. Thus the condition (27) is not fulfilled. 
Hence to stabilize the system, Az is not determined by its functional relation 
to Aq, but left as a free parameter and determined in the same way as Aq- 
Even for Ah, where its final value is exactly known, it is more advantageous 
to let it as a free parameter to stabilize the numerical algorithm. The starting 
values for A G , A z , and A H are chosen such that the infrared solution and 
the intermediate momentum solution fit continuously after the initial c^- have 
been set. After each global Newton run, the coefficients are refitted. Thus the 
Newton iteration is termed a micro-cycle while the fix-point iteration of the 
coefficients is termed macro-cycle. 

This approach is feasible, as it turns out that the radius of convergence for 
the coefficients is significantly larger than for the Chebychef coefficients. Aq 
is strongly constrained by the vanishing of the inverse ghost dressing function 
in the infrared, thus determining Az and Ah by their exact solutions. In other 
cases possibly another global iteration for the coefficients would be necessary. 
Note that during the iteration the values of the Ap can change by several 
orders of magnitude, even if initial and final value are close by. Furthermore 
as (27) is not fulfilled with the initial values it makes no sense to run the 
global Newton method until it converges. Therefore after a number of steps 
Nn the Newton method is interrupted, if it did not yet converge, to make a 
new macro-cycle step and then start a new micro-cycle. After several macro- 
cycle steps the Newton method already converges, i.e. \E\ < 5q for a fixed 5q, 
in less than N N steps. In the present case it was possible to choose N N fixed. 
In more instable solutions an adaptive choice may be appropiate. 

In addition, still slow convergence is encountered. To deal with it leads to the 
introduction of super-cycles: it turns out that convergence is extremely slow if 
e is chosen already quite small from the beginning. In this case the presented 
initial values are far from the final ones, and only very small steps can be made 
[23]. Since the infrared regime and the finite momenta regime do separate in 
the system under investigation, it is possible to choose an e larger than the final 
matching scale. The solutions then are not continuous differentiable at e and 



Special thanks to Claus Feuchter and Burghard Griiter for a very inspiring dis- 
cussion on this topic. 



14 



the relation (27) is not fulfilled. Thus the independence of the coefficients Ap 
in the fitting procedure is important here as well. Using the previous solutions 
for a given e as the input for a smaller one, it was possible to approach the 
correct solutions stepwise with sufficiently good convergence. Indeed, regularly 
only the first steps of the micro-cycle in the first one or two macro-cycle 
steps utilized the global properties of the method and afterwards local and 
thus quadratic convergence in \E\ was achieved. This is repeated until (27) is 
satisfied to the required precision. 

It is possible to reduce the number of necessary super-cycle steps drastically 
in two ways. One is an improvement of the fit functions Dp and fp to include 
sub-leading corrections. This leads to the final form of equation (34) and the 
fits (35-37). The other is the introduction of damping in the determination of 
the A F . To this end the new values A' F providing a continuous fit between D F 
and }'fDp are computed, but instead setting the value of A" F in the next step 
to A' F , only a variable admixture of A' F to the original A F is chosen as 

*=*£^. (42) 
where d > is the damping parameter. 

To test whether the algorithm has succeeded, two different criterions for the 
micro- and macro-cycle are used. For a micro-cycle the absolute value of \E\ 
is used. The typical scale of terms is set by the tree- level term of order 1. 
Compared to this scale a mean deviation of the order 10~ 8 is achieved with 
the described numerical parameters. The largest deviation are usually found 
in the far infrared of equation (11), where diverging terms have to cancel and 
thus very high accuracy is necessary. 

If the available precision is not sufficient, the necessity for these cancellations 
is also a limiting factor for the algorithm. An example of this problem is 
encountered when ( 3 in case of the varying solution branch described 
in subsection 3.1. In this case, the divergencies of the two competing terms 
are softened, and thus very small values of e are necessary before only the 
infrared leading terms are relevant and thus the asymptotic infrared solution 
Dp is sufficient. If ( is sufficently close to 3, the required values of e are so 
small that the cancellations cannot be resolved with the available precision 
and the algorithm fails to converge. 

For the macro-cycle the largest rate of change of the A F is used as a measure 
of convergence. Typical when the rate of change is below a fixed value 5j of 
the order of a per mill the macro-cycle is completed. For the super-cycle the 
ratio AqA z / J((,eo) according to (27) is used 9 . Note that due to numerical 

9 The value of Ah is always so close to the final value that this is not a relevant 
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errors it is even with a machine precision of 10 16 only possible to fulfill (27) 
to the order of 10 -5 . 



3.4 Notes on implementation and optimization 

While the implementation of the macro- and super-cycles is straightforward, 
there are some subtleties concerning the implementation of the micro-cycle. 
Once the values of the A F are fixed in the macro-cycle, calculating the cur- 
rent value of the right hand sides of the DSEs is straightforward numerically. 
However, as the evaluation of the exponentials in the fits of the F is numeri- 
cally costly, it is much more efficient to calculate all equations for a given pj 
simultaneously, as here the same set of Gauss-Legendre points can be used. 

This issue is much more important when performing the Newton step. To 
determine the descent direction for the coefficients Cy, it is necessary to calcu- 
late the Jacobian. As this Jacobian contains the derivatives of the equations 
(39) with respect to the coefficients 10 c^-, each element contains essentially a 
two-dimensional integration. Calculating therefore each element one-by-one is 
extremely inefficient. It is therefore advantageous to evaluate all elements with 
fixed pj simultaneously for all derivatives, corresponding to a complete column 
in parallel. This parallelism is implemented by performing one addition of the 
Gauss-Legendre integration for each of the equations simultaneously by sum- 
ming in an appropriate vector. Therefore it is also not possible to use standard 
Gauss-Legendre integrators from libraries, as they usually do not allow in a 
two dimensional array efficient parallel summation for a vector of integrals 
with appropriate treatment of input functions like the dressing functions. 

The performance can be enhanced further by close inspection of the system 
under consideration. In the case of equations (9-11) it can be exploited that 
inside the integrals one of the dressing function does not depend on the in- 
tegration angle. Excluding it from the angular integration also improves the 
cancellation of divergencies, thus reducing the required number of angular in- 
tegration points. The same applies to factors of Tj in the calculation of the 
Jacobian matrix. 

Compared to the time necessary to calculate the Jacobian, all other elements 
are completely subleading. Therefore it is no advantage to use an approximate 
inversion routine to speed up the calculation of the inverse Jacobian for the 
Newton step. The gain in speed does not outweigh the drawback in precision. 
It is furthermore useful to store the inverted Jacobian and not to recalculate 

figure of merit. 

10 Note that the derivative of exp(cijTj) w.r.t. the Cij is Tj exp(cyTj). Hence no 
costly numerical derivatives of the dressing functions are necessary. 
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Table 2 

The physical parameters in the example run. 
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Table 3 

The numerical parameters in the example run. 



Initial Aq 


1.06894 


Final A G 


0.838007 


Initial Az 


0.99109 


Final A z 


13.4111 


Initial Ah 


1.01411 


Final A H 


1.01404 


Initial ratio A G A Z / J((, ec) 


0.12035 


Final ratio A 2 G A Z / J((, e G ) 


1.0019 


Initial ratio A^/m^ 


1 


Final ratio A^/m 2 


0.99993 


Initial \E\ 


3.6 x 10 5 


Final \E\ 


2.5 x 10~ 9 


Final rate of change of Aq 


7.9 x 10~ 6 


Number of micro-cycles 


25 



Table 4 

The initial and final values for the infrared coefficients together with some run 
statistics. The run time was about 20 minutes on a 2.8 GHz Intel Xeon processor. 

it during the backtracking, as it does not change. Using all these optimization 
gives a factor of roughly 30-60 in speed. 



3.5 Example calculation 



Here the results for an example calculation will be presented. The physical 
parameters are listed in table 2. The numerical parameters can be found in 
table 3 and the initial values along with the final values and run parameters 
in table 4. The final values for the coefficients Qj are not displayed. Initially 
they have been set to 0. Using these parameters results have been obtained in 
a single super-cycle step. 

The results are presented in figures 1-3 for the ghost, gluon, and Higgs dressing 
functions, respectively. Along the various final fitting functions in (34) are 
shown. Together with the results in table 4 this demonstrates the abilities of 
the method in a most direct way. Note especially that the mismatch measured 
by \E\ is decreased by 14 orders of magnitude. Up to 18 orders have also been 
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Fig. 2. The gluon dressing function compared to the infrared, intermediate, and 
ultraviolet fitting functions. 

achieved. Further examples can be found in refs. [7,8,9]. 



4 4-dimensional Yang-Mills theory 



In case of the renormalized system (15-17), the main qualitative change is 
the appearance of the renormalization constants, defined in equations (18-24). 
Especially the infrared properties of the po = modes do not change [7,9]. In 
the ultraviolet the coefficients in (28) change, but this effect can be ignored. 
In fact the subleading term is completely dropped and only the leading term 
is (and must be) retained. The po ^ modes become constant both in the 
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Fig. 3. The Higgs dressing function compared to the infrared, intermediate, and 
ultraviolet fitting functions. 

infrared and the ultraviolet. 



If the renormalization constants would only appear explicitly as in equations 
(15-17), then the most direct implementation would be to treat them on equal 
footing with the remaining integrals. Although this would add 3N ch integrals 
when evaluating the Jacobian, this would be a subleading effect. However, 
the constants cannot be kept the same during a complete Newton iteration. 
Although the system is finite, the effect of renormalization is still large and 
the would-be divergencies would lead to solutions which are incorrect, thus 
slowing the convergence significantly if not leading to divergence. 

In addition, this is not the only effect. Due to the truncations made [7,9], it is 
necessary in nearly all integrals to replace the factor F(q ,q)F(p Q ± q ,p± q) 
by" 

( F( *- « - zkj) K ± * - zd**>) ' (43) 

Therefore the renormalization constants appear in a non-trivial combination, 
which are especially cumbersome when derivatives with respect to the Cij are 
performed. Although this problem could also be dealt with by determining 
the renormalization constants and their derivatives once for each N ch set of 
equations, a more direct approach is possible. For a single Newton step, it is 
viable to keep the renormalization constants fixed. Then afterwards new renor- 
malization constants are determined and the next Newton step is performed. 
Therefore each Newton iteration is split into sub-iterations, each of one step 
length, which is repeated in a micro-cycle for N N steps or until convergence 
has been achieved. Due to the global convergence of the Newton method this 

11 Z3 is again the generic wave-function renormalization constant of the dressing 
function F. 
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is permissible and quite efficient, as the renormalization constants not change 
drastically within one Newton step. This strategy may not be possible in cases 
where these constants do change appreciably. In the present case, this is the 
most efficient solution. This is important, as now 3NN ch equations have to 
be solved. With this method, a system of iV = 20 and N ch = 60 (yielding a 
Jacobian with about 10 7 entries) converges without further prior knowledge 
within a few days on a 2.8 GHz Intel Xeon. 



5 Conclusions 



In this work a globally convergent method was presented to solve coupled 
sets of nonlinear integral equations. The example investigated are DSEs for 
Yang-Mills theories. Equipped with qualitative knowledge of the infrared and 
ultraviolet properties only, it was possible to solve the system without further 
prior knowledge with an acceptable efficiency. Most crucial to this efficiency 
is the careful optimization of the calculation of the Jacobian inside the global 
Newton method, which is at the heart of the algorithm. 

The approach showed global convergence, thus extending earlier local algo- 
rithms, and therefore provides the abilities to solve such sets of equations 
without referring to other methods to provide good initial guesses. It will 
therefore be useful in the further investigations of DSEs and similar mathe- 
matical problems. 

However, the method relies crucially on knowledge of the asymptotic behavior. 
Attempts to at least remove the need for the infrared behavior have failed so 
far, although many different methods have been tried, ranging from expansion 
in suitable characteristic functions up to genetic algorithms. In many physical 
cases these asymptotics are known, but especially in gauge theories this is 
usually not yet the case in general gauges, see e.g. [24]. Therefore it remains 
desirable to construct such a numerical continuum method. 
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A Integral kernels and tadpoles of the 3-dimensional System 



The integral kernels in (9-11) are obtained using standard methods [1,7,8]. 
They depend only on the magnitude k and q and the relative angle 6 between 
their two momentum arguments. The kernel in the ghost equation (9) is given 
by 

, , q 2 sin 3 (fl) 

T[ ,q) ~ (k 2 + q 2 -2kqcos6) 2 ' 

The contributions in the Higgs equation (10) are 



,r n \ 2q 2 sm 3 (9) t , 

N ^ q) = - { k 2 + q 2 + 2kqLe ? (A ' 1} 
, T „ \ 2sin 3 (#) 

N ^ = - k * + q * + 2k qmS » - < A ' 2) 
The kernels in the gluon equation (11) are finally 



R(k,q) 



((C ~ ljfcg cos(fl) - q 2 + (q 2 cos 2 (9)) sin 9 
2k 2 (k 2 + q 2 + 2kq cos I) 
((C- l){k 2 + Akq cos 6) -4g 2 + 4g 2 Ccos 2 (#))sin# 



Mi(A; ' g) " Ak 2 {k 2 + q 2 + 2kqcos9) (A ' 3) 

m ^ *> = 4 ^ + r+2w + ^ - ^ - ^ 

+8(C - 3) (A; 2 + q 2 )kq cos 9 + (8(q 4 + (C + 7)k 4 
+4(5C - l)A; 2 g 2 ) cos 2 (fl) + 4(4(g 2 + (C + 3)A; 2 ) cos 3 (fl) 

+4(k 2 q 2 cos 4 ^)). (A.4) 

The modified gluon vertex is introduced by multiplying M T with 

(A(q, G, Z)A(q + k, G, Z)A(k, G, Z)) 5 ** (A.5) 
where A(q, G, Z) is given by 

A(q,G,Z) = l —— — . (A.6) 

Z(q)G(q) 2+ ^ 

The tadpoles used read: 
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TGG = -^f 2 J dqd9(R D (p, q )(G( q )G(p + q )-A 2 G q- 2 ^(p + qr 2 ^) 



+R 3 (p,q)(G(q)G(p + q)-Alq- 2 ^(p + q)~ 2e o - ^ 
+M TD (p,q)Z(q)Z(p + q)). (A.7) 



Here the kernel R has been split into a convergent, ^-independent part Rq, its 
divergent part Rd, and a remaining part R3 as 

R = R + (C - 3)# 3 + Rd- (A.8) 

The divergent parts are in general isolated by 

R D (p,q) = lim R{p,q). (A.9) 

M T £) contains the divergent part of M T , including the modifications due to 
the dressed 3-gluon- vertex (A. 5). The other tadpole in the gluon equation is 
given by 



J dqde(M LD (p,q)H(q)H(p + q)y (A.IO) 



Again, Mld contains the divergent part of Ml- In the Higgs equation, the 
tadpoles are set to 

iffG , rpHH _ 9 '3C A mh 

p 2 4tc 

A detailed account for the construction of the tadpoles is given in ref. [7]. 



rpHG _|_ rpHH _ "I'll ^ ^ 
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